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Low  Frequency  Propagation 
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This  report  examines  the  utility  and  limitations  of  the  integral- 
equation  representation  of  ground-wave  propagation  over  nonuniform  terrain. 
Emphasis  is  on  frequencies  between  20  kHz  and  200  kHz.  The  one-dimensional 

\ version  of  the  Integral  groundwave  equation  is  subject  to  errors  caused  by: 

44-  topographic  irregularities  near  the  great-circle  propagation  path; 
f2)  finite  ground  conductivity;^-)-  nonuniformities  in  the  earth's  electrical 
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properties ;A4y  an  approximate  Integration  to  reduce  the  dimension  of  the 
equation  from  two  to  one.  Each  of  these  errors  is  quantified,  and  the 
types  of  terrain  to  which  the  Integral  ground-wave  equation  is  applicable 
are  defined. 

A method  of  numerical  solution  is  developed  and  used  to  obtain  results 
for  the  special  case  of  a smooth,  uniform,  spherical  earth.  These  results 
are  compared  In  detail  with  numerical  results  obtained  from  the  widely 
used  residue-series  representation  of  ground -wave  propagation.  The  agree- 
ment between  the  two  methods  is  shown  to  be  excellent.  Graphical  results 
are  given  for  the  ground-wave  attenuation  function. 


jlCuairr  CL amifiCation  Of  r»*i*  O 


1.  I have  reviewed  the  Final  Report  by  Pacific  Sierra  Research 
Corporation,  titled  PROPAGATION  OF  THE  LOW  FREQUENCY  GRCUNDWAVE 
OVER  NONUNIFORM  TERRAIN.  The  report  describes  work  performed 
under  contract  to  RADC/EEP,  and  represents  the  product  of  that 
work. 


2.  The  report  describes  the  analysis  of  errors  associated 
with  various  formulations  which  mathematically  describe  the 
propagation  of  low  frequency  waves  over  the  earth's  surface. 

Each  of  the  errors  is  quantified,  and  the  types  of  terrain 
to  which  the  integral  ground-wave  equation  is  applicable  are 
defined.  A method  of  numerical  solution  is  developed  and  used 
to  obtain  results  for  the  special  case  of  a smooth,  uniform, 
spherical  earth.  The  results  of  this  study  apply  directly  to 
in-house  efforts  being  conducted  by  the  Propagation  Branch  (EEP) 
to  support  the  ESD  LORAN  SPO's  requirement  for  an  extremely 
accurate  means  of  predicting  LF  groundwave  propagation  parameters 


3.  All  aspects  of  the  work  contracted  for  by  RADC/EEP  have  been 
completed  and  described  in  the  subject  report  by  Pac.  Sierra 
Res.  Corp.,  and  meets  the  approval  of  this  office.  The  efforts 
by  Pacific  Sierra  have  resulted  in  information  of  importance  to 
our  understanding  of  the  propagation  of  long  radio  waves  over 
non-uniform  terrains. 


Paul  A.  Kossey,\EEP/X' 
Contract  Monitor 
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SUMMARY 


This  report  examines  the  utility  and  limitations  of  the  integral- 
equation  representation  of  ground-wave  propagation  over  nonuniform 
terrain.  Emphasis  is  on  frequencies  between  20  kHz  and  200  kHz.  The 
one-dimensional  version  of  the  integral  ground-wave  equation  is  subject 
to  errors  caused  by:  1)  topographic  irregularities  near  the  great-circle 
propagation  path;  2)  finite  ground  conductivity;  3)  nonuniformities  in 
the  earth's  electrical  properties;  4)  an  approximate  integration  to 
reduce  the  dimension  of  the  equation  from  two  to  one.  Each  of  these 
errors  is  quantified,  and  the  types  of  terrain  to  which  the  integral 
ground-wave  equation  is  applicable  are  defined. 

A method  of  numerical  solution  is  developed  and  used  to  obtain 
results  for  the  special  case  of  a smooth,  uniform,  spherical  earth. 

These  results  are  compared  in  detail  with  numerical  results  obtained 
from  the  widely  used  residue-series  representation  of  ground-wave 
propagation.  The  agreement  between  the  two  methods  is  shown  to  be 
excellent.  Graphical  results  are  given  for  the  ground -wave  attenuation 
function. 
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I.  INTRODP  CT ION 


The  propagation  velocity  of  low-frequency  groundvaves  Is  subject  to 
perturbations  from  nonuniformities  In  either  the  topographic  or  electrical 
properties  of  the  terrain.  A sufficiently  accurate  theory  of  propaga- 
tion over  irregular  terrain  would,  in  principle,  make  it  possible  to 
correct  position  errors  that  such  velocity  perturbations  cause  on  low- 
frequency  radio  navigation  systems. 

Historically,  two  theoretical  treatments  of  groundwave  propagation 
have  evolved:  the  residue  series  of  Van  der  Pol,  Bremer,  Horton,  and 
Fock  (see  e.g.,  Bremer,  1958);  and  the  integral  equation  approach  (see 
e.g.,  Bufford.,  1952 ; and  Feinberg , 1959).  The  residue  series,  being 
more  amenable  to  analytic  solution,  has  been  the  foundation  of  most  pre- 
vious results.  However,  although  useful  for  analysis  of  propagation 

* 

over  a uniform — or  a piecewise  uniform  — earth,  the  residue  series  is 
awkward  for  analysis  of  propagation  over  continuously  varying  terrain. 

For  the  latter  conditions,  numerical  solution  of  the  Integral  ground- 
wave  equation  appears  the  most  fruitful  approach. 

The  Integral  equation  is  based  on  impedance  boundary  conditions, 
which  are  approximate.  Therefore,  regardless  of  the  numerical  accuracy 
of  its  solution,  the  classical  version  of  the  integral  equation  cannot 
provide  accuracy  better  than  that  inherent  in  the  impedance  boundary 
condition.  Beyond  implementation  of  these  boundary  conditions,  a number 
of  additional  approximations  are  usually  made,  with  the  result  chat  the 
computationally  simplest  versions  of  Che  relevant  equations  are  subject 
to  the  most  stringent  conditions  on  the  types  of  terrain  to  which  they 
are  applicable. 

Accordingly,  a main  purpose  of  this  report  is  to  quantify  the  accuracy 
of  several  forms  of  the  integral  ground-wave  equation,  thereby  ascertain- 
ing the  types  of  terrain  to  which  they  may  be  used  to  within  specified 
error  tolerances.  Attention  is  restricted  to  frequencies  between  20  kHz 
and  200  Hz,  especially  to  the  LORAN-D  frequency  at  100  kHz.  The  approach 

*An  example  of  a piecewise  uniform  earth  would  be  two  or  more  uniform 
regions  separated  by  an  abrupt  boundary,  such  as  shoreline. 
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taken  Is  to  derive  expressions  for  errors  due  to  the  approximate  treat- 
ment of  1)  finite  earth  conductivity,  2)  terrain  curvature,  and  3)  non- 
uniformities in  electrical  properties.  These  error  terms  arise  in  two 
pieces  in  the  derivation  of  the  one-dimensional  form  of  the  integral 
equation:  1)  use  of  impedance  boundary  conditions;  2)  use  of  the  method 
of  stationary  phases,  or  something  nearly  equivalent,  to  perform  an 
integration  over  the  coordinate  transverse  to  the  great-circle  propaga- 
tion path,  thereby  reducing  a two-dimensional  equation  to  a one-dimensional 
equation. 

Use  of  Impedance  boundary  conditions  is  essential  to  the  derivation 
of  the  classical  integral  equation,  and  the  resulting  inaccuracies  must 
be  considered  inherent  to  the  formulation.  The  reduction  from  a two- 
dimensional  to  a one-dimensional  integral  equation,  however,  is  a 
simplification  that  could  be  forgone  at  the  expense  of  an  order-of- 
magnltude  Increase  in  difficulty  in  obtaining  numerical  solutions.  To 
determine  when  such  an  increase  in  difficulty  would  Indeed  provide  a 
commensurate  Increase  in  overall  accuracy,  we  compare  the  error  terms 
due  to  the  approximate  transverse  integration  with  the  ones  due  to  the 
impedance  boundary  conditions. 

Section  II  gives  the  relevant  versions  of  the  integral  ground-wave 
equation,  and  rank-orders  the  various  error  terms,  which  are  derived  in 
Appendices  A and  B.  Section  III  gives  numerical  results  comparing  two 
versions  of  the  one-dimensional  Integral  equation  with  each  other  and 
with  the  residue  series.  Section  IV  presents  conclusions;  Appendix  C 
derives  the  Integral  equation  in  polar  coordinates  for  a smooth,  round 
earth;  and  Appendix  D outlines  the  procedures  used  to  obtain  numerical 
solutions. 
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II.  INTEGRAL  EQUATIONS  FOR  GROUNDWAVE  PROPAGATION 

We  begin  by  outlining  the  steps  required  to  derive  the  classical 
one-dimensional  integral  equation  of  Hufford  (.1952) . Although  the  final 
• result  is  simply  the  well-known  Eq.  (11)  of  Huf ford's  original  paper, 

the  intermediate  steps  reveal  somewhat  more  general  forms.  Moreover, 
we  Identify  the  points  at  which  critical  approximations  are  made,  and 
quantify  the  accuracy  of  these  approximations.  Finally,  we  give  a form 
of  the  integral  equation  that  is  somewhat  more  accurate  than  Hufford 's 
for  the  special  case  of  a smooth,  round  earth. 

HUFFORD* S EQUATION 

Although  awkward  for  the  special  conditions  of  a smooth,  round  earth , 
rectangular  coordinates  (x,y,z)  are  the  most  convenient  when  the  shape 
of  the  earth's  surface  cannot  be  given  a simple  analytic  form.  We  assume 
1)  that  the  transmitter  is  located  at  the  origin,  and  let  ?(x,y)  denote 
the  deviation  of  the  earth's  surface  from  the  plane  z«0;  and  2)  that 
the  receiver  is  located  in  the  vertical  plane,  y*0.  Below,  we  also  use 
an  integration  point,  Q,  which  is  on  the  earth's  surface  and  has  coordina- 
tes x,y,c(x,y),  and  a receiver  point,  P,  with  coordinates  Xq,?(x^).  Other 
quantities  used  below  are  r^,  r^,  and  r^,  which  are  straight-line  dis- 
tances between  the  origin  and  P,  the  origin  and  Q,  and  Q and  P,  respect 
tively.  (Appendix  B gives  expressions  for  rQ,  r^,  and  r^.) 

The  refractive  index  of  the  earth,  n^,  is  given  by 

n^  - < + io/we  , (1) 

S 0 

. where  k Is  the  dielectric  constant  of  the  earth,  o is  the  conductivity  of 

the  earth,  m is  the  angular  frequency  of  the  wave,  and  is  the  vacuum 
dielectric  permitlvlty.  All  parameters  in  Eq.  (1)  can  be  spatially  non- 
uniform,  although  the  validity  of  the  forthcoming  equations  depends  on 
these  nonuniformities  falling  within  constraints  given  below. 

For  certain  smooth,  symmetric  surfaces,  (e.g.,  planes  or  spheres), 
the  Hertz  potential  for  a vertical  electric  dipole  has  only  a sing. a 


4 


component  oriented  normal  to  the  surface.  In  such  Instances,  only  three 

(E  ,E  ,H  ) field  components  are  excited  that  can  be  calculated  from  a 
x z y 

single  potential  function,  <j».  For  arbitrary  rough  surfaces,  the  Hertz 
vector  Is  not  oriented  In  the  normal  direction,  and  all  six  field  com- 
ponents are  excited.  To  avoid  this  complexity,  the  derivation  of  the 
Integral  ground-wave  equation  Is  based  on  the  assumption  that  the  surface 
Is  "so  smooth"  that  the  Hertz  vector  Is  composed  of  essentially  a single 
component  oriented  In  the  z direction. 

Quantification  of  the  error  involved  in  this  approximation  is  dif- 
ficult. Intuitively,  one  would  expect  it  to  be  valid  provided  the  in- 
clinations of  the  earth's  surface  with  respect  to  its  average  level  are 
small.  More  rigorously,  if  y is  the  angle  between  the  z axis  and  the 
normal  to  the  earth's  surface,  the  error  will  be  roughly  the  amount  by 
which 

cosy  • 


differs  from  unity. 

Subject  to  the  error  given  by  Eq.  (2),  <i»  satisfies  the  wave  equation 

(72  + 1c2)*  - t (3)' 

where  t is  the  source  function,  k is  the  free-space  wave  number,  and  a 
time  dependence  e has  been  assumed. 

The  second  major  approximation  is  use  of  Impedance  boundary  con- 
ditions, which  can  be  stated  in  the  form 


fi  - -«c«*  , 

where  n is  the  upward  normal  to  the  earth's  surface  and 


5 


5 = 1/n 


provided  Chat  n Is  large.  The  accuracy  of  Che  impedance  boundary  con- 


dition  (discussed  in  detail  in  Appendix  A)  is  summarized  below.  Use  of 
the  impedance  boundary  conditions  is  essential  to  the  derivation.  It 
permits  9i|>/3n  to  be  expressed  in  terms  of  4*  on  the  surface,  which,  in 


turn,  permits  use  of  Green's  Theorem  to  convert  Eq.  (3)  to  an  Integral 

* 

equation  involving  integration  over  the  earth  s surface.  The  procedure. 


described  by  Hufford,  gives  the  result 


*<"'  - V?)  + Hr/’W  ^ [i+(1+cj)Tr  • <« 


where  i|*n  is  the  potential  that  would  dxist  if  the  earth  were  not  present. 


By  letting 


♦Q( Q)  - Const. 


and  defining  W by 


<KQ)  - 2W(QU  (Q)  , 


we  find 


. (7) 


More  accurate,  albeit  much  more  complicated,  versions  of  Eq.  (4)  can 
be  derived  from  the  results  of  Rytov  {1940) . 
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In  Eq.  (7),  U is  an  attenuation  function  that  accounts  for  the  fact  that 
the  earth  is  not  flat  and  does  not  have  infinite  conductivity.  Note 
that  H ■ 1 if  S’  ■ 0 (o  - <■)  and  3r2/3n  - 0. 

Equation  (7)  is  the  most  general  form  of  the  Integral  equation,  being 
subject  only  to  the  limitations  of  the  impedance  boundary  conditions  and 
the  assumption  of  gentle  departures  from  a plane  earth.  Being  a two- 
dimensional  integral  equation,  however,  Eq.  (7)  is  quite  expensive  to 
solve  numerically.  A major  simplification  (carried  out  in  Appendix  B) 
transforms — subject  to  some  restrictions — Eq.  (7)  into  a one-dimensional 
integral  equation.  Formally,  this  transformation  involves  using  the 
method  of  stationary  phases  to  perform  analytically  the  integration  over 
the  coordinate  transverse  to  the  propagation  path.  Physically,  this 
transformation  implies  that  only  regions  within  the  first  Fresnel  zone 
significantly  affect  the  received  signal.  The  resulting  equation  is 


W(Xq)  - 1-e 


W(x)(6  + 3r2/3n)e 


Ut(rl+r2’r0) 


(8) 


which  is  the  classic  form  derived  by  Hufford  (1952). 

ACCURACY  OF  HUFFORD'S  EQUATION 

Equation  (8)  is  an  order-of -magnitude  simpler  to  solve  than  Eq.  (7), 
but  is  less  accurate  because  of  errors  incurred  in  the  approximate  trans- 
verse integration.  It  is  therefore  Important  to  quantify  the  accuracy 
of  these  equations  to  determine  whether  Eq.  (7)  is  sufficiently  more 
accurate  (or  more  general)  than  Eq.  (8)  to  warrant  the  considerable 
additional  computational  complexity.  Moreover,  it  is  importfut  to  esta- 
blish the  limitations  on  Eqs.  (7)  and  (8),  thereby  determini  q the  types 
of  terrain  to  which  each  may  be  applied  to  achieve  some  specified  accuracy. 
Accordingly,  Table  1 summarizes  the  first-order  correction  terms  to  each 
of  the  main  approximations.  These  correction  terms  (derived  in  Appendices 
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A and  B)  denote  the  order  of  magnitude  of  the  errors*  involved  in  each 
approximation.  In  Table  lt  denotes  the  local  radius  of  curvature  of 
the  boundary. 


Table  1 

ACCURACY  OF  MAIN  APPROXIMATIONS 


Approximation 

Fractional  Error 

I.  Hertz  vector  normal  to  surface 

0C/Sx)2  + (H/3y)2 

II.  Impedance  boundary  conditions 

a.  Finite  conductivity 

“eo/ff 

b.  Surface  curvature 

[*Vt]'1/2/R0 

c.  Nonuniform  conductivity 

win[m 

III.  Stationary  phase  integration 

a.  Stationary  point  at  ya  0 

<*/»7>J|y.0 

b.  Asymptotic  series 

l/kx0 

Table  1 shows  that  errors  due  to  use  of  1)  impedance  boundary  con- 
ditions for  finitely  conducting  madia  (Ila)  and  of  2)  the  asymptotic 
expansion  of  the  stationary  phase  integration  (Illb)  are  the  most  funda- 
mental in  the  sense  that  they  are  nonzero  even  for  a plane,  uniform  earth. 

The  other  error  terms  depend  on  the  degree  of  terrain  nonuniformity. 

A simple  conclusion  regarding  the  relative  accuracy  of  approximations 
Ila  and  Illb  cannot  be  made,  because  one  depends  on  ground  conductivity, 
whereas  the  other  depends  on  the  length  of  the  propagation  path.  Table 
Al  (p.  30)  gives  numerical  values  for  the  term  Ila,  and  shows  that  better 
than  1-percent  accuracy  is  obtained  at  LF/VLF  provided  that  a > 10  mhos/m. 

* 

Roughly  speaking,  the  percentage  error  associated  with  each  approxima- 
tion can  be  estimated  by  multiplying  the  fractional  errors  of  Table  1 by  100. 
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Poor  accuracy  is  obtained  for  Greenland  ice,  i.e.,  where  o a 10  ^ mhos/m. 

If  the  stationary-phase  error-term  Illb  is  less  than  the  impedance  error- 
term  Ila,  no  degradation  in  accuracy  (for  a plane  earth)  is  caused  by 
the  transformation  of  Eq.  (7)  to  Eq.  (8),  and  any  additional  accuracy 
achieved  by  solving  the  two-dimensional  equation  would  be  spurious.  Con- 
versely, if  the  term  Illb  exceeds  Ila,  additional  accuracy  is  obtained 
by  dealing  with  the  complexity  of  the  two-dimensional  equation.  Comparison 
of  these  terms  shows  that  the  stationary  phase  integration  causer  no  sub- 
stantial degradation  in  accuracy  provided  that  x^  exceeds  a characteristic 
distance,  1,  given  by 


(9) 


Table  2 gives  1 for  various  conductivities  and  a frequency  of  100  kHz. 

_3 

These  results  show  that  if  the  conductivity  is  10  mhos/m  or  less,  the 

one-dimensional  integral  equation  is  essentially  as  accurate  as  the  two- 

dimensional  equation,  provided  that  the  pathlength  exceeds  about  100  km. 

_2 

For  a conductivity  of  10  tnhos/m,  the  two-dimensional  equation  is  more 
accurate  than  the  one-dimensional  one,  unless  the  pathlength  exceeds 
about  900  km.  For  seawater  (a  * 4 mhos/m) , the  two-dimensional  equation 
is  far  more  accurate  than  the  one-dimensional  one  for  all  realistic  path- 
lengths  . 


Table  2 

VALUES  OF  l AT  100  kHz  FOR  SEVERAL  CONDUCTIVITIES 


a (mhos/m) 

4 

H* 

O 

1 

KJ 

IQ'3 

•n 

1 

o 

•H 

A (km) 

3.5  x 105 

8.6  x 102 

90 

0.9 
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Before  evaluating  the  other  expressions  in  Table  1,  note  that  their 
derivation  involves  power  series  expansions,  and  that  these  expansions 
are  valid  only  when  their  magnitude  is  less  than  unity.  Caution  must 
also  be  exercised  in  interpreting  the  error- term  I for  a spherical  earth. 

* For  a smooth,  round  earth,  it  is  easy  to  show  that 

2 

(3p/3x)2  + (3p/3y)2  - , (10) 

l-(x/a)Z 

...  - 

which  correctly  indicates  that  Eq.  (8)  becomes  very  inaccurate  as  the  pro- 
pagation pathlength,  x,  approaches  an  earth  radius,  a.  This  unnecessary 

I inaccuracy,  however,  is  due  to  Hufford's  treatment  of  the  earth's  curvature 

as  a perturbation  to  a plane  earth.  The  situation  is  remedied  by  deriving 
the  integral  equation  in  spherical  coordinates  (Appendix  C and  below) , 
which  causes  the  appropriate  Hertz  vector  to  be  rigorously  normal  to  the 
surface  for  a smooth,  round  earth.  Thus,  the  proper  Interpretation 
* of  such  error  terms  as  1 and  Ilia  should  treat  C as  the  departure  of  the 

terrain  contour  from  the  average  surface  contour  of  the  earth;  l.e.,  C 
should  Include  hills,  etc.,  but  not  the  earth's  curvature,  which  can  be 
accurately  accounted  for. 

For  a frequency  of  100  kHz,  Table  A2  (p.  31)  shows  that  errors  due  to 
surface  curvature  (expression  Ila  in  Table  1)  are  about  an  order-of- 
magnltude  greater  than  those  due  to  finite  conductivity  (expression  lib), 

-3 

for  conductivities  of  10  mhos/m  or  more,  and  • 1 la.  For  R^  ■ 100  m, 
the  error  terms  due  to  curvature  effects  exceed  10  percent  for  normal 
ground  conductivities.  Nonthelsss,  in  this  case,  it  is  better  to  account 
for  hills  via  Eqs.  (7)  or  (8)  than  to  leave  them  out  of  the  analysis 
entirely.  Note  that  by  setting  RQ  ■ a in  expression  lib,  it  follows  that 
the  Impedance  error  caused  by  normal  earth  curvature  is  extremely  small. 

For  Greenland  ice,  the  results  of  Appendix  A show  that  errors  due 
to  finite  conductivity  effects  at  100  kHz  are  so  large  that  discussion 
of  other  error  sources  is  academic.  Ve  do  not  have  adequate  data  to 
evaluate  the  term  IIIc,  which  accounts  for  nonuniformities  in  conductivity. 
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It  Is  easily  shown,  by  evaluating  expressions  I or  Ilia,  that  errors 
caused  by  terrain  gradients  should  be  less  than  1 percent  for  grades  of 
5 percent  or  less,  and  less  than  10  percent  for  grades  of  15  percent  or 
less.  Much  steeper  grades  would  essentially  totally  destroy  the  accuracy 
of  either  Eq.  (7)  or  Eq.  (8). 

In  summary,  for  very  smooth  terrain  where  the  error  terms  (Table  1) 

I,  lib,  c,  and  Ilia  are  much  smaller  than  the  terms  Ila  and  Illb,  the 

two-dimensional  integral  equation  is  much  more  accurate  than  the  one- 

_2 

dimensional  equation  only  for  ground  conductivities  of  10  mhos/m  or 
more.  For  much  lower  ground  conductivities,  the  accuracy  of  the  Impedance 
boundary  conditions  is  sufficiently!  poor  that  nothing  additional  Is 
really  lost  by  resorting  to  the  approximate,  one-dlmenslonal  Integral 
equation.  Once  grades  of  5 percent  or  more,  or  terrain  features  with 
radius  of  curvature  of  1 km  or  less  are  encountered,  the  one-dlmen- 
slonal equation  might  as  well  be  used,  since  the  accuracy  of  the 
stationary  phase  integration  is  no  worse  that  that  of  the  other 
approximations . 


INTEGRAL  E 


ROUND  EARTH 


Equation  (8)  can  be  used  to  calculate  W for  the  case  of  a smooch, 
round  earth  provided  chat  the  propagation  path  does  not  exceed  a mega- 
meter or  so.  As  shown  by  Eq.  (10),  the  errors  caused  by  departure  of 
the  earth's  surface  from  the  plane  z * 0 can  be  substantial  for  longer 
propagation  paths.  Such  errors  can  be  avoided  by  rederiving  the  Integral 
equation  In  spherical  coordinates,  and  using  Che  radial  component — rather 
than  the  z-component — of  the  Hertz  potential.  In  this  Instance,  the 
error  terms  I and  Ilia  In  Table  1 vanish,  whereas  the  term  lib  Is  extre- 
mely small  for  ■ a.  The  result  Is  that  an  Integral  equation  can  be 
obtained  that  Is  essentially  as  accurate  for  a uniform  spherical  earth 
as  is  Eq.  (8)  for  a uniform  plane  earth.  In  addition,  this  "spherical" 
Integral  equation  provides  a consistent  basis  of  comparison  with  results 
calculated  using  the  residue  series,  which  rigorously  accounts  for  a 
spherical  earth. 

* 

The  resulting  one-dlmenslonal  Integral  equation  Is 

* 

See  Appendix  C for  detailed  derivation. 
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«<.„> 


,tB  IT 

sin  s/2a 


sin  s/a 

80-s 

sin  Sq/ a sin 


[s  + .in  ^-]  . 


aQ-«1  2ika jsin  s/ 2a -sin  Sg/2a-**sin  - j 


where  s and  s^  denote  great-circle  distances  on  the  earth's  surface. 

By  using  rQ  - 2a  sin  Sq /2a,  etc.,  in  the  exponent  of  Eq.  (11),  and 
noting  that 


3r2/3n  - sinjj~-J  , 


it  follows  that  Eqs.  (8)  and  (11)  agree  to  the  extent  that  the  approximation 


«0  *0 
sin  — = — 


Is  valid;  l.e.,  the  fractional  disagreement  between  Eq.  (8)  and  Eq.  (11) 

2 

is  of  order  (Sg/a)  , which  arises  because  Eq.  (8)  is  subject  to  errors 
of  the  order  of  magnitude  indicated  by  Eq.  (10),  whereas  Eq.  (11)  is  not. 
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III.  NUMERICAL  RESULTS  FOR  A UNIFORM.  SPHERICAL  EARTH 

Equation  (11)  (p.  11)  is  the  most  accurate  form  of  the  one-dimensional 
integral  equation  for  a smooth,  spherical  earth,  and  is  therefore  solved 
numerically  to  obtain  the  results  given  in  this  section.  Specifically, 
the  attenuation  function,  W,  is  computed  as  a function  of  distance  for 
frequencies  between  20  kHz  and  200  kHz,  and  conductivities  between 
4 mhos/m  and  2x10  ^ mhos/m.  (Appendix  D outlines  the  numerical  methods 
used.) 

Although  somewhat  less  accurate  than  Eq.  (11)  for  a uniform,  spherical 
earth,  Huf ford's  Integral  equation  (Eq.  (8))  is  convenient  for  analyzing 
propagation  over  irregular  terrain.  Accordingly,  as  a partial  check  on 
relative  accuracy,  we  also  solve  Eq.  (8)  (p.  6)  numerically  and  compare 
the  results  with  those  obtained  from  Eq.  (11).  In  solving  Eq.  (8),  we 
used  the  full  expressions  for  rQ,  r^,  r^,  and  3r2/3n  (e.g.,  rQ-2a  sin  sQ/ 
2a,  etc.)  rather  than  the  expression  in  powers  of  s/a  as  x/a  used  in 
Huf ford's  (1952)  example. 

As  an  accuracy  check  on  both  Eqs.  (8)  and  (11),  detailed  comparisons 
with  results  given  for  the  residue  series  by  Wait  and  Howe  (1956)  are 
made.  Care  must  be  exercised  in  making  these  comparisons,  because  we 
have  defined  the  attenuation  function  by  (see  Eq.  (C-9)) 


ikr„ 


!ji(s0)  ©<  W(sQ)e 


(12) 


whereas,  a different  attenuation  function,  W,  is  defined  by 


iks. 


*(8Q) 


W(s0)e 


(13) 


in  the  residue  series.  Thus,  since  the  Hertz  potential,  ,j, , must  be  the 
same  in  both  treatments,  the  phase  of  W given  by  Wait  and  Howe  must  be 
corrected  by  a factor 


k(s 


0 


k 


- 2a  sin 


» 


(14) 
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before  comparison  with  our  results.  This  adjustment  In  the  residue  series 
results  has  been  made  In  the  comparisons  given  below. 

To  make  the  comparisons  as  quantitative  as  possible,  we  use  a digital 
rather  than  a graphical  format.  Tables  3 through  16  give  the  results; 
the  number  of  significant  figures  given  corresponds  to  the  numerical 
accuracy  that  we  used  In  solving  Eqs.  (8)  and  (11).  The  results  labeled 
"spherical  int.  eq."  correspond  to  Eq.  (11);  those  labeled  "Hufford  Int. 
eq."  correspond  to  Eq.  (8);  those  labeled  "residue  series"  are  taken 
from  Ualt  and  Howe,  adjusted  according  to  Eq.  (14).  Following  Halt  and 
Howe,  we  used  an  effective  earth  radius  of  4a/3  (e.g.,  8S00  tan)  to  account 
for  atmospheric  refraction.  The  sensitivity  of  the  results  to  the  choice 
of  effective  radius,  which  Is  crude  at  low  frequencies.  Is  examined  below. 
For  purposes  of  comparison  with  Walt  and  Howe,  we  used  k ■ 0 In  the  cal- 
culations, but  do  not  advocate  doing  so  In  general.  Also,  the  rather 
unusual  distances  at  which  the  results  are  given  were  chosen  so  that 
comparison  with  Walt  and  Howe  could  be  made. 

We  discuss  first  the  results  at  100  kHz,  since  this  frequency  Is 
of  more  practical  Interest  than  the  others.  Tables  3 through  6 give 
these  results,  and  show  that  both  Che  Hufford  equation  and  Eq.  (11)  agree 
with  the  residue-series  results  to  within  one  tenth  of  a degree  of  phase 
for  distances  out  to  600  kilometers.  To  put  this  accuracy  In  context, 
note  that,  at  100  kHz,  one  tenth  of  a degree  of  phase  corresponds  to  a 
distance  of  less  than  a meter.  Moreover,  even  this  minute  disagreement 
is  due  to  roundoff,  and  would  have  been  smaller  had  more  significant 
figures  been  presented  In  the  tables. 

At  distances  of  1200  km  or  more,  small — but  noticeable— differences 

appear  between  the  results  of  Eq.  (8)  and  Eq.  (11)  and  the  residue  series. 

Equation  (11)  agrees  somewhat  more  closely  with  the  residue  series  than 

does  Eq.  (8).  This  behavior  is  to  be  expected,  because  Eq.  (8)  is 

2 

accurate  only  to  order  (.Sq/a)  , which  becomes  appreciable  (0.08  at  2400 
km)  at  the  larger  distances.  We  have  no  way  of  knowing  whether  the 
disagreement  between  Eq.(ll)  and  the  residue  series  at  2420  km  (Table  3 ) 

* 

Numerical  accuracy  pertains  to  the  precision  of  the  methods  used 
to  solve  the  equations,  and  has  nothing  whatever  to  do  with  the  accuracy 
of  the  equations  themselves,  which  Is  discussed  In  Sec.  II  and  Appendices 
A and  B. 


Table  3 


■I  ■ #»■  n 


W AS  COMPUTED  FROM  EQ.  (11)  FOR  f * 100  kHz  AND  a = 4 mhos/m 


Distance  (km)  60.6 


Phase  U (deg) 
Residue  series 


Phase  W (deg) 
Spherical  Int  eq 


Difference  0 


121 


4.3 


4.3 


0 


242 


10.9 


10.8 


0.1 


606 


47.8 


47.8 


0 


1211 


199.5 


199.2 


0.3 


2420 


1148.2 


1151.8 


(3.6) 


Amplitude  W 

Residue  series 

0.983 

0.952 

0.869 

0.576 

Amplitude  W 

Spherical  Int  eq 

0.982 

0.951 

0.868 

0.575 

Difference 

O^Ol 

0.001 

0.001 

6.001 

0.223 


0.223 


0 


0.024 


0.023 


0.001 


Table  4 

W AS  COMPUTED  FROM  EQ/  (8)  FOR  f =■  100  kHz  AND  a = 4 mhos/m 


Distance  (km) 

60.6 

121 

242 

606 

1211 

Phase  U (deg) 

Residue  series 

2.0 

4.3 

10.9 

47.8 

199.5 

Phase  W (deg) 

Hufford  Int  eq 

2.0 

4.3 

10.9 

47.9 

198.8 

Difference 

0 

0 

0 

(0.1) 

Amplitude  W 

Residue  series 

0.983 

0.952 

0.869 

0.576 

0.223 

Amplitude  W 

Hufford  Int  eq- 

0.982 

0.951 

0.868 

0.575 

0.222 

Difference 

0.001 

0.001 

0.001 

0.001 

0.001 

IS 

Table  5 


W AS  COMPUTED  FROM  EQ.  (11)  FOR  f = 100  kHz  AND  a = 10'2  mhos/m 


Distance  (km) 

60.6 

121 

242 

606 

1211 

Phase  W (deg) 
Residue  series 

20.1 

30.1 

47.7 

109.2 

297.0 

Phase  VI  (deg) 
Spherical  Int  eq 

20.1 

30.0 

47.6 

109.2 

296.9 

Difference 

0 

0.1 

0.1 

0 

0.1 

Amplitude  U 

Residue  series 

0.969 

0.927 

0.828 

0.531 

0.206 

Amplitude  VI 
Spherical  Int  eq 

0.969 

0.927 

0.828 

0.531 

0.206 

Difference 

0 

0 

0 

0 

0 

Table  6 

W AS  COMPUTED  FROM  EQ.  (8)  FOR  f 3 100  kHz  AND  o - 10-2  mhos/m 


Distance  (km)  60.6  121  242  606  1211 


Phase  W (deg) 
Residue  series 

20.1 

30.1 

• 47.7 

109.2 

297.0 

Phase  W (deg) 
Hufford  int  eq 

20.1 

30.0 

47.7 

109.3 

296.9 

Amplitude  W 

Residue  series 

0.969 

0.927 

0.828 

0.531 

0.206 

Amplitude  W 

Hufford  int  eq 

0.968 

0.926 

0.829 

0.531 

0.206 

Difference 

0.001 

0.001 

(0.001) 

0 

0 

16 

Table  7 

U AS  COMPUTED  FROM  EQ.  (11)  FOR  f - 20  kHz  AND  a * 4 mhos/m 


Distance  (km) 

60.6 

121 

242 

606 

1211 

Phase  W (deg) 
Residue  series 

■a 

1.6 

4.2 

17.5 

61.6 

Phase  W (deg) 
Spherical  int  eq. 

829 

4.1 

17.4 

61.5 

D1 f ference 

0 

0.1 

0.1 

0.1 

0.1 

i I 

Amplitude  W 

Residue  series 

0.992 

0.978 

0.939 

0.779 

0.497 

Amplitude  W 
Spherical  Int  eq 

0.991 

0.977 

0.937 

0.777 

0.495 

D1 f ference 

0.001 

0.001 

0.002 

0.002 

0.002 

Table  8 

W AS  COMPUTED  FROM  EQ.  (8)  FOR  f » 20  kHz  AND  o = 4 mhos/m 


Distance  (km) 

60.6 

121 

242 

606 

1211 

Phase  U (deg) 
Residue  series 

■a 

1.6 

4.2 

17. ‘ ’ 

61.6 

Phase  W (deg) 
Hufford  Int  eq 

■a 

1.6 

4.2 

17.5 

61.9 

Difference 

0 

0 

0 

0 

(0.3) 

i 

Amplitude  W 

Residue  series 

0.992 

0.978 

0.939 

0.779 

0.497 

Amplitude  W 

Hufford  Int  eq 

0.992 

0.977 

0.938 

0.778 

0.495 

Difference 

0 

0.001 

0.001 

0.001 

0.002 

Table  9 


U AS  COMPUTED  FROM  EQ.  (11)  FOR  f 3 20  khz  AND  a = 10‘‘  mhos/m 


Distance  (km) 


Phase  W (deg) 
Residue  series 


Phase  VI  (deg) 
Spherical  int  eq 


Difference 


Amplitude  W 

Residue  series 

0.993 

0.978 

0.938 

0.780 

0.504 

Amplitude  W 
Spherical  int  eq 

0.991 

0.976 

0.936 

0.778 

0.502 

D1 f ference 

0.002 

0.002 

0.002 

0.002 

0.002 

Table  10 

W AS  COMPUTED  FROM  EQ.  (11)  FOR  f * 20  kHz  AND  a * 10"3  mhos/m 


Distance  (km) 


Phase  W (deg) 
Residue  series 


Phase  W (deg) 
Spherical  int  eq 


Difference 


Amplitude  W 
Residue  series 


1211 


119.8 


119.5 


.3 


0.987  0.967  0.920  0.750  0.477 


?p"er1cealWint  .q  °'987  °-567  °-s2°  ».«7 


Difference 


0.001 


* 1 ^ 
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Table  11 

W AS  COMPUTED  FROM  EQ.  (11)  FOR  f =*  50  kHz  AND  a » 4 mhos/m 


Distance  (km) 

60.6 

121 

242 

606 

1211 

Phase  W (deg) 
Residue  series 

1.2 

2.8 

7.1 

30.5 

117.5 

Phase  W (deg) 
Spherical  Int  eq 

1.2 

2.8 

7.1 

30.5 

117.4 

Difference 

0 

0 

0 

0 

0.1 

I 

Amplitude  W 

Residue  series 

0.988 

0.966 

0.905 

0.675 

0.338 

Amplitude  U 
Spherical  Int  eq 

0.987 

0.965 

0.904 

0.674 

0.337 

Di f ference 

0.001 

0.001 

0.001 

0.001 

0.001 

Table  12 

W AS  COMPUTED  FROM  EQ.  (11)  FOR  f - 50  kHz  AND  a =*  10"2  mhos/m 


Distance  (km) 

60.6 

121 

Hi 

606 

1211 

Phase  W (deg) 
Residue  series 

10.3 

. 15.6 

25.0 

60.8 

164.4 

Phase  W (deg) 
Spherical  int  eq 

10.3 

15.6 

25.0 

60.8 

164.2 

Difference 

0 

0 

0 

0 

0.2 

I 

Amplitude  W 

Residue  series 

0.985 

0.960 

0.896 

0.667 

0.342 

Amplitude  W 
Spherical  int  eq 

0.985 

0.960 

0.896 

0.667 

0.342 

Difference 

0 

0 

0 

0 

0 
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Table  13 

W AS  COMPUTED  FROM  EQ.  (11)  FOR  f * 50  kHz  AND  a » 10-3  mhos/m 


Distance  (km) 

60.6 

121 

242 

606 

1211 

Phase  W (deg) 
Residue  series 

30.8 

44.6 

66.4 

126.5 

263.5 

Phase  W (deg) 
Spherical  Int  eq 

30.8 

44.6 

66.4 

126.5 

263.5 

Difference 

0 

0 

0 

0 

0 

i 

Amplitude  W 

Residue  series 

n 

0.899 

0.79C 

0.505 

0.213 

Amplitude  W 
Spherical  int  eq 

0.957 

0.903 

0.795 

0.509 

0.211 

Difference 

(0.005) 

(0.004) 

(0.005) 

(0.004) 

B 

Table  14 

W AS  COMPUTED  FROM  EQ.  (11)  FOR  f - 200  kHz  AND  a » 4 mhos/m 


Distance  (km)  60.6  121  242  606  1211 


Phase  W (deg) 
Residue  series 

Phase  W (deg) 
Spherical  Int  eq 


Difference 


Amplitude  W 

Residue  series 

0.976 

0.932 

0.820 

0.462 

0.130 

Amplitude  W 
Spherical  Int  eq 

0.975 

0.931 

0.819 

0.461 

0.129 

Difference 

0.001 

0.001 

0.001 

0.001 

0.001 

Table  15 

W AS  COMPUTED  FROM  EQ.  (8)  FOR  f - 200  kHz  AND  o » 4 mhos/m 


Distance  (km) 


Phase  U (deg) 
Residue  series 


Phase  VI  (deg) 
Hufford  Int  eq 


Difference 


606 

1211 

77.3 

350.9 

77.5 

353.9 

(0.2) 

(3.0) 

Amplitude  W 

Residue  series 

0.976 

0.932 

0.820 

0.462 

Amplitude  W 

Hufford  Int  eq 

0.975 

0.931 

0.819 

0.461 

Difference 

0.001- 

0.001 

0.001 

0.001 

Table  16 

W AS  COMPUTED  FROM  EQ.  (11)  FOR  f » 200  kHz  AND  a - 10'2  mhos/m 


Distance  (km) 


Phase  W (deg) 
Residue  series 


Phase  W (deg) 
Spherical  Int  eq 


Difference 


199.0  546.7 


199.0  544.7 


Amplitude  W 

Residue  series 

0.921 

0.834 

0.668 

Amplitude  U 
Spherical  Int  eq 

0.922 

0.836 

0.670 

Difference 

(0.001) 

(0.002) 

(0.002) 

0.002 
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Is  due  Co  Che  approximations  made  in  Chls  report,  or  numerical  impreci- 

slons  In  Wale  and  Howe's  resulcs.  Also,  Wale  and  Howe  did  noc  give  Che 

exacC  value  ChaC  Chey  used  for  Che  earch's  radius.  However,  even  for 

Che  worse  case  shown  (a  ■ 4 mhos/m,  s ■ 2420  km),  Che  dlsagreemenc  In 

o 

phase  Is  3.6°,  which  corresponds  Co  a dlscance  of  only  30  mecers. 

The  remaining  cables  (7  ehrough  16)  furcher  confirm  Che  general 
conclusions  drawn  above.  For  dlscances  up  Co  600  km,  Che  agreemenc 
among  Eq.  (8),  Eq.  (11),  and  Che  residue  series  Is  vlrCually  exacC.  For 
greacer  dlscances,  Che  agreemenc  Is  sclll  excellenC  buC,  as  expecCed, 

Eq.  (11)  agrees  sllghcly  more  closely  wich  Che  residue-series  resulcs 
Chan  does  Eq.  (8). 

The  resulcs  of  Tables  3 Co  16  are  sufflclenCly  close  Co  chose  of 

WalC  and  Howe  ChaC  a graphical  presenCaClon  here  would  add  nochlng  new. 

WalC  and  Howe,  however,  did  noc  presenc  resulcs  for  propagaclon  over  ice, 

nor  were  Chey  able  Co  obCain  sacisfaccory  convergence  of  che  residue 
-3 

series  for  o ■ 10  mhos/m  and  frequencies  of  100  kHz  and  200  kHz.  For 
Chese  condiCions,  Cherefore,  we  presenc  graphical  resulcs  (Figs.  1 
ehrough  4) . 

Figure  1 gives  Che  ampllCude  and  phase  of  W for  a frequency  of 

_3 

100  kHz  and  a ground  conducCiviCy  of  10  mhos/m.  Resulcs  are  shown 
for  boch  Che  "4/3"  earch  used  by  WalC  and  Howe,  and  a "normal  earch"  of 
radius  6372  km.  The  resulcs  for  Chese  Cwo  assumed  earch  radii  agree 
quiCe  closely,  alchough  noeiceable  differences  do  occur  aC  ranges  of 
several  hundreds  of  kilomecers.  For  example,  aC  a range  of  600  km, 

Che  selecCion  of  effecCive  earch  radius  can  Influence  Che  calculaced 
phase  by  more  chan  20  degrees,  which  corresponds  Co  a posicion  uncercalncy 
of  abouc  170  mecers.  Since  che  use  of  an  effecCive  earch  radius  co 
accounc  for  aCmospherlc  refracclon  is  crude,  chls  20-degree  difference 
beCween  Che  "4/3"  and  "normal"  earchs  muse  be  regarded  as  a sore  of 
uncercalncy,  which  far  exceeds  che  machemacical  uncertainties  associated 
wich  Eqs.  (8)  and  (11).  In  ocher  words,  Che  accuracy  of  the  equations 
seems  Co  be  far  better  Chan  chls  input  Co  Che  equations. 

Figure  2 gives  che  amplitude  and  phase  of  W versus  distance  for  a 

_3 

frequency  of  200  kHz  and  c ■ 10  mhos/m;  and  Figs.  3 and  4,  che  amplitude 
and  phase  of  W for  various  VLF/LF  frequencies,  and  electrical  properties 


normal  Earth 


Amplitude  and  phase  of  W for  f - 100  kHz  and  o - 10~J  mhoa/m 


r 


(Bap)  AA  9»Md 

ili§l*iS§§88S8o_ 


propagation  over  Greenland  ice 


■ 
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corresponding  Co  nominal  Greenland  ice  (e.g.,  < ■ 6,  and  a » 2 x 10  ^ 
mhos/m) . Given  Che  poor  accuracy  of  Che  impedance  boundary  condiCions 
for  low  conducCiviCy  (Table  1A) , Che  precision  of  Che  resulCs  in  Figs.  3 
and  4 is  noC  high.  As  expecced.  Figs.  1 Chrough  4 show  chac  Che  aCCenua- 
cion  and  phase  shifc  increase  as  Che  frequency  increases  and  as  Che 
ground  conducCiviCy  decreases. 


! 
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IV.  CONCLUSIONS 

For  a smooch,  uniform  earch,  Hufford's  one-dimensional  integral 
equation  gives  essentially  exact  agreement  with  the  results  of  the  resi- 
due series  for  frequencies  between  20  kHz  and  200  kHz,  and  propagation 
distances  up  to  several  hundreds  of  kilometers.  For  propagation  dis- 
tances greater  than  1000  km  to  1500  km,  the  accuracy  of  Hufford's  equation 
degrades  somewhat,  and  a modified  one-dimensional  equation — expressed  in 
polar  coordinates — provides  slightly  better  agreement  with  the  residue 
series. 

For  a nonuniform  earth  having  terrain  undulations,  the  two-dimension- 
al integral  equation  exhibits  errors  due  to:  1)  assuming  that  the  Hertz 
vector  is  essentially  normal  to  the  surface;  2)  use  of  Impedance  boundary 
conditions.  The  one-dimensional  integral  equation  incurs  additional 
errors  due  to  the  approximate  ^valuation  of  an  integral  over  the  coordi- 
nate transverse  to  the  propagation  path. 

_2 

For  ground  conductivities  greater  than  about  10  mhos/m,  or  for 
pathlengths  less  than  about  100  km,  the  two-dimensional  Integral  equation 
is  inherently  much  more  accurate  than  the  one-dimensional  one,  provided 
that  the  earth  is  fairly  smooth.  However,  this  additional  accuracy  could 
be  unnecessary,  because — as  was  the  case  for  the  perfectly  smooth, 

• spherical  earth — the  accuracy  of  the  one-dimensional  equation  could  be 
adequate. 

_3 

For  ground  conductivities  less  than  about  10  mhos/m,  or  for  rela- 
tively rough  terrain,  the  Inherent  accuracy  of  the  two-dimensional 
equation  is  really  no  better  than  that  of  the  simpler  one-dimensional 
version.  This  behavior  occurs  because  errors  due  to  the  assumption  of 
Impedance  boundary  conditions  and  a normally  oriented  Hertz  vector  are 
at  least  as  large  as  those  due  to  the  approximations  made  in  reducing  the 
two-dimensional  equation  to  the  classical  one-dlmenslonal  form.  More 
specifically,  for  these  unfavorable  terrain  characteristics,  the  accuracy 
of  the  two-dimensional  equation  is  degraded  to  the  extent  that  no 
additional  penalty  is  paid  for  performing  the  approximate  transverse 
integration. 
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Appendix  A 

ACCURACY  OF  IMPEDANCE  BOUNDARY  CONDITIONS 

Use  of  Impedance  boundary  conditions  (Eq.  (4),  p.  4)  Is  essential  to 
the  derivation  of  the  Integral  equation  for  the  groundwave  attenuation 
function.  Therefore,  the  accuracy  of  even  an  exact  solution  of  the  full- 
fledged  two-dimensional  Integral  equation  (Eq.  (7),  p.  5)  is  no  better 
than  the  accuracy  of  the  impedance  boundary  conditions.  The  applicability 
of  these  boundary  conditions  has  received  detailed  attention  by  numerous 
authors  (e.g.,  Rytov , 1940 ; Leontoviah , 1944;  Brermer,  1958;  Feinberg* 
1959;  Senior , 1961),  and  the  details  of  their  treatments  need  not  be 
repeated  here.  However,  to  quantify  the  inherent  accuracy  of  Eq.  (8) 
for  the  frequencies  and  terrain  of  interest  here,  this  appendix  briefly 
summarizes  the  formulas  for  the  correction  terms  to  the  impedance 
approximation . 

Impedance  boundary  conditions  are  accurate  for  highly  conducting, 
uniform  media  having  flat  boundaries.  Errors  are  thus  incurred  if  1)  the 
medium  is  imperfectly  conducting,  2)  the  electrical  properties  of  the 
medium  are  spatially  nonuniform,  and  3)  the  boundary  of  the  medium  is  not 
flat.  We  consider  each  of  these  errors  below. 

ERRORS  DUE  TO  FINITE  CONDUCTIVITY 

Validity  of  impedance  boundary  conditions  rests  on  the  fact  that  for 
a highly  conducting  earth,  the  refracted  wave  is,  according  to  Snell's 
law,  propagated  in  a direction  nearly  normal  to  the  earth's  surface. 

Even  for  a plane  earth,  deviation  of  the  wave  normal  in  the  earth  from 
the  normal  to  the  surface  causes  computational  errors.  Rytov  (.1940) 
and  Leontovlch  (1944)  showed  that,  for  a vertically  polarized  wave  and 
a uniform  plane  earth,  the  fractional  error  incurred  by  use  of  impedance 
boundary  conditions  is  of  order 


^"Fractional  error"  is  defined  as  the  ratio  of  the  neglected  terms  to 
the  retained  term  in  a series  representation. 
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where  we  have  assumed  chat  o/u>e  » k,  where  < Is  the  dielectric  constant 

0 

of  the  ground.  For  situations  where  o/u£q  is  not  much  greater  than  k, 

the  error  is  of  order  1/k.  Since  k is  only  of  order  10  or  less  for  most 

types  of  ground,  Eq.  (A-l)  applies  for  all  situations  where  the  accuracy 

of  Impedance  boundary  conditions  is  better  than  about  10  percent. 

To  quantify  the  error  term  (A-l) , we  give  its  numerical  values  for 

conductivities  and  frequencies  of  Interest  in  Table  Al.  These  values 

show  that,  for  most  cases,  use  of  impedance  boundary  conditions  is  valid 

in  the  sense  that  the  fractional  error  incurred  is  much  less  than  unity; 

l.e.,  the  percent  error  is  no  greater  than  about  10  percent.  If  very 

_2 

high  accuracy — say,  1 percent  (fractional  error  of  10  ) or  better — is 

required,  however,  the  impedance  boundary  conditions  are  inadequate  in 
the  20  kHz  to  200  kHz  range  for  Greenland  ice,  which  has  a conductivity 
of  about  10  3 mhos/m. 


Table  Al 

VALDES  OF  uCq/o  FOR  VARIOUS  FREQUENCIES  AND  CONDUCTIVITIES 


^\^o(mhos  /m) 
f (kHz) 

4 

O 

1 

10'3 

10“5 

20 

2.8  x 10~7 

1.2  x 10"4 

1.1  x 10~3 

1.1  x 10'1 

50  j 

6.9  x IQ-7 

2.9  x 10'4 

2.9  x 10~3 

2.8  x 10-1 

100 

1.4  x 10*6 

5.8  x 10-4 

5.6  x 10"3 

5.6  x 10'1 

200 

2.9  x 10-6 

1.1  x 10-3 

1.2  x 10-2 

1.2 

ERRORS  DUE  TO  SURFACE  CURVATURE 

The  errors  caused  by  curvature  of  the  boundary  have  been  calculated 
by  Rytov  (1940)  and  Senior  (1961).  Physically,  the  validity  condition 
is  that  the  local  radius  of  curvature,  R^,  of  the  surface  be  large  com- 
pared with  the  skindepth  of  the  wave  in  the  earth.  Mathematically,  the 
expression  for  the  fractional  error  Incurred  by  using  impedance  boundary 
conditions  in  the  presence  of  a curved  boundary  is 
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(*V£>'1/2/r0 


(A-2) 


To  quantify  this  error  term.  Table  A2  gives  values  for  the  expression 

3 

(A-2)  for  various  conductivities  and  frequencies  and  Rg  • 10  m.  One  way 
to  interpret  Table  A2  is  that  the  values  shown  are  the  fractional  errors 
due  to  hills  (for  example)  having  a radius  of  curvature  of  1 km.  The 
error  caused  by  hills  with  a 100-m  radius  of  curvature  would  be  10  times 
as  large  as  shown  in  Table  A2,  whereas  that  due  to  hills  with  a 10-km 
curvature  would  be  one  tenth  as  large. 

Table  A2 

VALUES  OF  (iru„af)"1/2/Ra  FOR  VARIOUS  CONDUCTIVITIES 
0 0 

AND  FREQUENCIES  AND  RQ  - 103m 


cr(mhos/m) 


f (kHz) 


1.8  x 10 


1.1  x 10 


8.0  x 10 


5.6  x 10 


3.6  x 10 


2.2  x 10 


1.6  x 10 


1.1  x 10 


1.1  x 10 


7.1  x 10 


5.0  x 10 


3.6  x 10 


7.1  x 10 


5.0  x 10 


3.6  x 10 


Table  A2  also  shows  that,  with  regard  to  curvature  effects,  the 
accuracy  of  the  impedance  boundary  condition  degrades  as  the  frequency 
decreases.  This  behavior  is  different  than  that  shown  for  finite-con- 
ductivity effects  (Table  Al) , where  the  accuracy  degrades  as  the  fre- 
quency is  increased.  Thus,  if  the  frequency  is  too  high,  the  Impedance 
approximation  fails  because  the  ground  refractive  index  is  too  small  to 
refract  the  wave  into  a nearly  normal  direction;  whereas  if  the  fre- 
quency is  too  low,  the  approximation  falls  because  the  sklndepth  can 
become  comparable  with  the  radius  of  curvature  of  the  surface.  Comparison 
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of  Tables  A1  and  A2  shows  Chat,  except  for  very  smooth  ground,  curvature 
effects  Induce  larger  computational  errors  than  finite-conductivity 
effects. 

The  errors  shown  In  Table  A2  apply  where  the  entire  propagation  path 
Is  characterized  by  undulations  having  characteristic  dimensions  of  a 
kilometer,  and  would  be  smaller  If  only  part  of  the  path  contained  such 
irregularities.  Also,  so  long  as  the  fractional  errors  shown  In  Table  A2 
are  less  thin  a few  tenths.  It  is  better  to  account  for  surface  undula- 
tions using  the  Impedance  method,  than  to  neglect  them  entirely.  On 

this  basis,  we  estimate  that  for  a frequency  of  100  kHz,  and  the  average 
-3  -2 

ground  (o  - 10  to  10  mhos/m) , the  Impedance  method  of  treating  curved 
terrain  should  be  useful,  provided  that  the  characteristic  dimensions  of 
the  undulations  are  at  least  100  meters.  For  smaller  values  of  Rq,  the 
error  terms  become  comparable  with  the  retained  terms. 

ERRORS  DUE  TO  NONUNIFORM  ELECTRICAL  PROPERTIES 

Rytov  (.1940) , Leontovich  (1944),  and  Senior  (1962)  have  evaluated  the 
errors  Incurred  by  using  Impedance  boundary  conditions  for  media  having 
nonuniform  electrical  properties.  Specifically,  for  k and  a that  are 
functions  of  x,  y,  and  z,  where  z is  the  vertical  coordinate,  they  showed 
that  the  fractional  error  Is  given  by 


1 

< 3z 


V (iTVlgfo)  l/f2 


(A-3) 


-1/2 

Since  (irUgfff)  is  simply  the  sklndepth  of  the  wave  In  the  ground,  the 
error  term  (A-3)  will  be  small  If  the  conductivity  undergoes  only  a small 
fractional  change  within  a sklndepth  of  the  surface.  The  expression  (A-3) 
is  obviously  extremely  small  for  seawater  because  the  sklndepth  Is  small 
and  the  medium  is  nearly  uniform.  At  a frequency  of  100  kHz,  we  see — by 
multiplying  the  values  In  Table  A2  by  1000 — that  the  sklndepth  In  normal 
ground  Is  several  tens  of  meters.  Thus,  for  Eq.  (A-3)  to  be  small,  the 
ground  must  be  nearly  homogeneous  to  a depth  of  tens  of  meters.  Precise 
evaluation  of  Eq.  (A-3)  must  use  data  on  shallow  conductivity-depth 
profiles  for  terrain  of  interest. 
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As  pointed  out  by  Senior  (.1962)  and  Leontovlch  (1944) , the  condition 
(A- 3)  depends  only  on  vertical  variations  In  conductivity,  even  through 
a was  assumed  to  have  lateral  Inhcmogeneitles  as  well.  Thus,  to  £lrst 
order,  nonuniformities  In  electrical  properties  Introduce  errors  only 
to  the  extent  that  the  wave  can  penetrate  to  a depth  where  the  con- 
ductivity departs  from  its  surface  value.  Errors  due  to  lateral  non- 
uniformities are  of  higher  order,  and  are  appreciable  only  where  the 
quantity 


(A-4) 


Is  very  large;  e.g.,  at  a coastline. 


I 
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Appendix  B 

VALIDITY  CRITERIA  FOR  ONE-DIMENSIONAL  INTEGRAL  EQUATION 


Appli tj.on  of  Impedance  boundary  conditions  and  Green's  theorem 

leads  to  the  following  type  of  Integral  equation  for  the  attenuation 
function,  W(P) 


W(P) 


- 1 + 


QW(Q)K(Q,P) 


(B-l) 


In  Eq.  (B-l),  the  integration  is  taken  over  the  surface.  A,  of  the  earth, 

Q is  an  integration  point  on  this  surface,  and  P is  the  receiver  point, 
which  we  assume  is  also  on  the  surface.  By  making  certain  approximations, 
the  two-dimensional  integration  over  the  surface.  A,  can  be  converted  into 
a one-dimensional  integration  along  the  terrain  between  transmitter  and 
receiver.  This  appendix  determines  the  accuracy  of  these  approximations. 

The  detailed  form  of  the  right-hand  side  of  Eq.  (B-l)  depends  on 
the  coordinate  system  used.  However,  the  validity  criteria  for  the 
approximate  integration  of  the  right-hand  side  of  Eq.  (B-l)  do  not  depend 
on  the  coordinate  system.  Therefore,  to  keep  the  algebra  as  simple  as 
possible,  we  establish  the  accuracy  of  the  approximate  integration  in 
rectangular  coordinates.  Appendix  C rederlves  the  integral  equation  in 
spherical  coordinates,  which  are  more  natural  to  propagation  over  a 
spherical  earth. 

Consider  a rectangular  coordinate  system  with  the  transmitter  at  the 
origin,  and  with  the  plane  defined  by  z • 0 oriented  perpendicular  to  the 
vertical  electric-dipole  transmitting  antenna.  We  let  the  deviation  of 
the  surface  of  the  earth  from  the  z*0  plane  be  given  by  C(x,y),  and 
assume  the  receiving  antenna  to  be  located  in  the  vertical  plane  defined 
by  y-0.  The  integration  point,  Q,  has  the  coordinates  x,  y,  C(x,y), 
and  the  following  relationships  hold: 


r 
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ro  “ x0  + 50  » co  “ c^xo'0^  » 


2 2 , 2 , 2.  . 

rx  - x + y + C (x,y)  , 


*2  " r*o’x^  + ^ + C?0_5(x,y)}2 


d2Q 


dxdy  ^1+fe)  + (ty) 


(B-2) 

(B-3) 

(B-4) 

(B-5) 


In  rectangular  coordinates,  Eq.  (B-l)  can  be  written  (see  Sec.  II  or 
Bufford,  1962) 


PcH 


r 2 ***( 

W(P)  - 1+1  p(Q)  e 

J 12 


(B-6) 


where 


F(Q)  -^W(Q)  [*+(l+^-)^] 


(B-7) 


By  using  the  relationships  (B-2)  through  (B-5),  Eq.  (B-6)  can  be' rewritten 


W(xQ)  - 1 + 


/-/ 


ikr 


dyf (x,y)e 


r-^-o 


(B-8) 


where 


•[*  +(i+^)£]  • 


(B-9) 


! 
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Thus,  Co  reduce  Eq.  (B-9)  Co  a one-dimensional  integral  equation,  our 
task  is  to  evaluate  the  integral 


I(x) 


m 

■f 


dyr(x,y)e 


ikr.h(x,y) 


(B-10) 


where 


rl  + r2 

h - — L ■ - 1 


(B-ll) 


For  transmission  paths  much  greater  than  X/2*  (about  0.5  km  at  100  kHz), 


krQ  » 1,  and  the  integral  (B-10)  is  of  the  classic  form  amenable  to  approxi- 
mate evaluation  by  the  method  of  stationary  phases  (e.g.,  Erdalyi,  19S6) . 

At  this  point  of  the  derivation,  other  authors  (e.g.,  Bufford , 1952)  have 
correctly  argued  that  the  integral  (B-10)  is  approximately  given  by  Che 
stationary-phase  formula,  which  is  well-known  and  can  be  written  down  by 
inspection.  This  formula  is,  in  fact,  the  leading  term  in  an  asymptotic 
series  representation  of  the  integral  (B-10) . To  quantify  the  accuracy  of 
this  term,  we  must  retain  and  evaluate  the  second-order  correction  terms. 

We  assume  that  h has  a stationary  point  at  y-y.,  given  by  the  equation 


h'  (y  ) - — ! 

'70J  3yl 


yyf 


- o 


(B-12) 


The  value  of  y^  will  be  found  below;  but,  for  now,  we  need  only  assume  that 
such  a stationary  point  exists.  Because  of  Eq.  (B-12),  the  power  series 
for  h becomes 


h(y)  - h+-j-(y-yQ)  +-g-(y-y0)  +—(y-yn)  + •••» 


(B-13) 


where  the  prime  denotes  differentiation  with  respect  to  y,  and  all  deri- 
vatives are  evaluated  at  ysyQ.  In  the  conventional  stationary  phase 
method,  only  the  first  two  terms  in  Eq.  (B-13)  are  retained — the  others 
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being  correction  terns  that  are  snail  If  krQ  » 1.  We  retain  these  higher- 
order  terns  and  assune — subject  to  a posteriori  justification— -that  they 
are  snail  enough  to  pernit  the  following  series  expansion  of  the  exponent 
In  Eq.  (B-10) : 

lkr0h  lkr0h(yQ)  lkxQh"  (y-yQ)2/2 


• 1+lkr, 


fh"*  / .3  . h""  . .4-1 1 

+_24-(y_y0)  Jj  * 


(B-14) 


By  similar  reasoning,  and  subject  to  similar  a posteriori  justification. 


we  write 


r r'(,,y0)(y-y0)-l 

r(i,7)  = r(x,y0)  |l+ r(  ) ■]  • 


(B-15) 


By  Inserting  Eqs.  (B-14)  and  (B-15)  Into  Eq.  (B-10),  and  noting  that  odd 
powers  of  (y-yg)  vanish  due  to  odd  symmetry,  the  Integral  becomes 


ikr  h(y  ) ? 

l(x)  :q  e 0 r(x,yQ)  Id  ye 


. [■,ritooh'"r' , ikto"'-i  .«] 

1 |_  6r  24  Jy“y05 


which  can  be  Immediately  Integrated  to  give 


. iri/4  ikrQh<y0 ) / 2 s \1/2rf  x 
Kx)  s e e r(X’y0) 

[L  lhmr' lh ,,n  ~[ 

2krQ(hM)2r  8kr0(h”)2J 


ikr0h"'  (y-yQ)  /2 


(B-16) 


(B-17) 
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where,  again,  a prime  denotes  y-di££erentiation  and  all  quantities  are 
evaluated  at  y-y^. 

Aside  £rom  the  question  o£  the  limits  an  the  x- Integration,  and  one 
or  two  minor  constraints,  insertion  o£  Eq.  (B-17)  into  Eqs.  (B-10)  and 
(B-8),  will  give  the  classical  one-dimensional  integral  equation  (e.g., 
Bufford,  1952 ; Bremer,  1958 ) provided  that  the  following  two  conditions 
are  satisfied! 


h(7Q)  a:  h(0) ; r^y^  s=  r^O)  etc.. 


r’ 


2kr0(h")2  F 


1 . 


8kr0(h")2 


« 1 


(B-18) 


(B-19a) 


(B-19b) 


Condition  (B-18)  aimply  states  that  tha  stationary  phase  point,  y^,  is  suf- 
ficiently close  to  the  plane  defined  by  ya0,  that  y^  s 0 may  be  sub- 
stituted in  all  relations.  This  condition  results  in  an  Integration  along 
the  line  between  transmitter  and  receiver.  Conditions  (B-19)  require 
that  itr^  be  large,  which  is  the  validity  requirement  for  the  stationary 
phase  integration.  To  quantify  the  accuracy  of  these  approximations,  we 
further  simplify  and  evaluate  the  correction  terms. 

By  Inserting  Eq.  (B-ll)  into  Eq.  (B-12),  performing  the  differentiation, 
and  using  a perturbation  expansion,  we  obtain  the  following  equation  for 
the  stationary  phase  point: 


(B-20) 

If  the  transverse  derivative  of  the  terrain  contour,  vanishes  at  y«0, 
then  7q"0,  and  we  recapture  the  classical  result  that  the  integrand  be 
evaluated  on  the  line  z»0,  y-0.  In  fact,  Eq.  (B-20)  shows  that  the 


_ [ri(C0-OT2c]C 


rl  + r2 
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stationary  point  is  displaced  from  the  line  z*0,  y ■ 0 by  an  amount  pro- 
portional to  the  lateral  gradient  of  the  terrain  contour,  C ' , evaluated 
on  the  line  between  transmitter  and  receiver. 

It  follows  from  Eq.  (B-17)  that  the  phase  of  I(x)  is  governed  by 


r0h(V  “ rl  + r2"r0  * <B_21) 

By  using  Eqs.  (B-2)  through  (B-4),  and  keeping  leading  terms  in  a power 
series  expansion  of  r^,  r^,  and  r^,  if  follows  that 


r0h(y0) 


2 2 2 2 2 

y0+r  yQ+(c0-c)  Cq 

2x  + 2(xq-x)  2^ 


(B-22) 


for  0 <,  x <_  Xq. 

From  Eq.  (B-20),  it  follows  that  the  order  of  magnitude  of  y^  is  eg', 
because  r^/^  + r^)  and  r^/fr^  + r^)  are  of  order  of  unity.  Therefore, 
the  fractional  phase  error  incurred  by  setting  yg  ■ 0 is 

0 [(£ 

To  the  same  order  of  accuracy,  the  y-derivative  term  in  Eq.  (B-5)  may  also 
be  neglected. 

Tedious,  but  straightforward  differentiation  of  h(y),  and  insertion 
of  ya0,  show  that  the  error  incurred  by  neglecting  the  terms  given  by 
(B-19)  is  of  order 


)21 
y - 0 / _ 


(B-23) 


h"w 

8kr0(h")2 


~ 0 (1/kTg) 

y • 0 


(B-24) 


i 


It  also  follows  from  direct  differentiation  that,  to  the  accuracy  indicated 
by  Eqs.  (B-23)  and  (B-24), 


r0h" (y-0) 


rl~l*  r2 
rlr2 


whence  Eq.  (B-17)  can  be  written 


(B-25) 


I(x) 


iri/4  lk(rl+r2_r0) 


(B— 26) 


where  all  quantities  are  evaluated  at  y-0. 

Comparison  of  Eq.  (B-26)  with  Eq.  (B-10)  shows  that  the  stationary- 
phase  integration  is  equivalent  to  multiplying  the  integrand  of  Eq.  (B-10) 
by  a lateral  distance.  Ay,  whose  magnitude  is  given  by 


Ay 


r 2ffrir2  i1/2  _ r xrir2  i1/2 

[lc(rl+r2)J  L(rl+r2iJ 


(B-27) 


which  is  essentially  the  width  of  the  first  Fresnel  tone.  Therefore,  the 
stationary-phase  integration  is  physically  equivalent  to  retaining  con- 
tributions from  transverse  distances  of  the  order  of  a Fresnel  zone. 
Satisfaction  of  condition  (B-18)  (or,  equivalently,  if  the  error  term 
(B-23)  is  small)  guarantees  that  the  terrain  will  not  vary  appreciably 
across  this  zone.  Combination  of  Eq.  (B-26)  with  Eqs.  (B-8)  through 
(B-ll)  gives  the  following  integral  equation  for  the  attenuation  function 


W(x0) 


1-  e 


-wi/4 


Wl 


dxyl  + Oc/3x)‘ 


ik(r,+r,-r.)  f 3r,l 
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Eq.  (B-22)  show*  that  the  exponent,  k(r1+r2-r0),  is  of  order  ItcVx  pro- 
vided that  0 <_  x <_  xQ.  Thus,  for  terrain  Irregularities  that  are  com- 
parable to  or  smaller  than  a wavelength,  the  exponent  In  Eq.  (B-24)  Is 


small  and  oscillates  slowly.  For  x < 0 or  x > xQ,  however,  the  exponent 
Is  of  order  2kx  or  21c(x-Xq),  which  oscillates  very  rapidly.  Thus,  to 
within  the  same  order  of  accuracy  as  the  stationary-phase  Integration 
(e.g.,  Eq.  (B-24)),  the  Integration  limits  In  Eq.  (B-28)  can  be  taken 
equal  to  0 and  Xq. 

Use  of  these  finite  limits  In  the  x-integratlon  converts  Eq.  (B-28) 
Into  virtually  Hufford's  classic  form.  However,  Hufford  made  two  addi- 
tional approximations,  which— although  consistent  with  those  described 

ebove  and  in  Sec.  II — are  not  really  necessary.  First,  Hufford  assumed 

2 2 
that  d Q ■ dxdy,  which  Is  tantamount  to  neglecting  (3g/3x)  In  Eq.  (B-28). 

Second,  he  assumed  that  terrain  Irregularities  were  sufficiently  gentle 

that,  except  In  the  exponent.  It  Is  permissible  to  set 


r0  * X0 


ri  « x 


r2  " X0-X 


(B-29) 


The  approximation  (B-29)  is  accurate  to  order  Cq/*q*  To  within  the  accuracy 
of  these  approximations,  Eq.  (B-24)  becomes  Eq.  (8)  (p.  6),  which  ts  identcal 
to  Eq.  (11)  In  Hufford's  (1952)  original  paper. 
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except  at  the  transmitter  location.  In  the  above  equations,  r Is  the 
distance  from  the  origin — in  this  Instance,  the  center  of  the  earth — and 


therefore  differs  qualitatively  from  rQ,  r^,  and  r2»  which  are  linear 


distances  between  transmitter  and  receiver,  transmitter  and  some  integra- 
tion point,  and  integration  point  to  receiver. 

Continuity  of  E and  H at  r ■ a implies  that 


«0T*I 


£tm,] 


e[r$] 


(C-5) 


Subject  to  the  validity  conditions  for  the  impedance  boundary  conditions 
given  in  Appendix  A, 


3/3r  [nji]  » . [r\p]  if  r < a 


and  the  boundary  condition  becomes 


yy  [riji]  • -ik6  [ rip  ] 


(C-6) 


which  is  simply  a spherical  coordinate  version  of  Eq.  (4)  (p.  4). 

Equations  (C-4)  and  (C-6)  may  be  used  in  conjunction  with  Green's 
theorem  to  obtain  a two-dimensional  integral  equation.  The  steps  are 
identical  to  those  used  by  Hufford  (1952)  for  the  plane  earth,  except 
that  here  the  volume  is  bounded  by  two  disconnected  surfaces — one  a large 
sphere  with  a radius  approaching  infinity,  the  other  coincident  with  the 
earth's  surface  except  for  an  infinitesimal  hemisphere  about  the  receiver 
location.  The  resulting  integral  equation  is 


*<r0)  “ 2¥q  + 


ikr. 


dAt|>(r^) 


(c- 


7) 
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where  Che  integral  is  taken  over  the  surface  of  the  earth  and  is  the 
free-space  Hertz  potential. 

We  let 


♦qC^)  - Const. 


(C-8) 


and 


♦ (r^)  - 2W(rx)  ^O^) 


(C-9) 


where,  as  above,  W denotes  an  attenuation  function  accounting  for  departures 
from  the  situation  where  the  earth  is  flat  and  of  Infinite  conductivity. 
Insertion  of  Eqs.  (C-8)  and  (C-9)  into  Eq.  (C-7)  gives  the  following  inte- 
gral equation  for  W. 


W(rQ) 


- 1 + 


r ik(r.+r,-r.) 
dAW(r.)  e 1 2 0 

x rlr2 


(c-10) 


Equation  (C-10)  is  formally  nearly  identical  to  Eqs.  (B-6)  and  (B-7) , 
the  extra  term  in  the  square  brackets  of  Eq.  (C-10)  resulting  from  the 
polar  spherical  coordinate  system.  For  a smooth,  round  earth,  however, 
explicit  functional  forms  can  be  given  for  dA,  r^»  r^,  etc.,  rather  than 
expressing  them  in  terms  of  the  unspecified  terrain  function,  ?,  and  its 
derivatives.  These  function  forms  permit  the  transverse  part  of  the 
surface  Integral  to  be  performed  without  having  to  use  a full-fledged 
stationary-phase  approximation. 

We  thus  write  (after  considerable  rearrangement) 

dA  - a2sinfl  d6  d*  , (C-ll) 


rQ  - 2a  sinflQ/2 


(C-12) 


9 


r^  • 2a  sine/ 2 


(C-13) 
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•'/la.  [1-CO80  cos6  - sin0  sin0n  cos$] 


1/2 


(C-14) 


3r 


1 1/2 
- — [l-cos0  cos0g-sin0  sin0Q  cos$] 

r - a /l 


(C-15) 


±.hi 

r2  3r 


x 

2a 


(C-16) 


r » a 


where  [9^,0]  and  [9,$]  are  the  angular  locations  of  the  receiver  and  the 
Integration  point  on  the  surface  of  the  earth.  Substituting  Eqs.  (C-ll) 
through  (C-15)  Into  Eq.  (C-10)  gives  the  following  form  for  the  integral 
equation 


W(0Q) 


sin0Q/2n  211ca[sln0/2  - sin0Q/2] 

sine/2  J e 


(C-17.) 


where 


I(k,0) 


i/2  ka[A - B cos<|>]1/2 


[A  - B cos$] 


1/2 


(C-18) 


and 


A ■ l-cos0  cos0q  , (C-19) 

B - sin0  »in0g  • (C-20) 

Mote  that  although  the  earth' 3 surface  is  azimuthally  symmetric,  r2 
depends  on  4 and  a transverse  Integration  (C-18)  must  be  performed. 
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Equation  (C-17)  is  much  simpler,  however,  than  a full-fledged  two-dimen- 
sional integral  equation  because  the  unknown  function,  W,  depends  only  on 
6.  The  transverse  integral,  I(k,8),  depends  only  on  known  quantities  and 
could  be  tabulated  via  numerical  Integration,  i.e.,  I(k,8)  is  a known 
function  of  8 that  need  not  be  determined  as  part  of  the  solution  process. 
Further,  note  that  the  only  approximation  used  to  derive  Eq.  (C-17)  is 
the  application  of  the  Impedance  boundary  conditions,  and  that  numerical 
solution  of  Eq.  (C-17) — in  conjunction  with  numerical  tabulation  of  l(k,8) — 
would  give  W(8)  to  the  accuracy  of  the  impedance  boundary  conditions. 
However,  as  is  evident  from  the  close  agreement  of  the  numerical  results 
given  in  Sec.  Ill  with  those  computed  by  other  authors  using  the  residue 
series,  an  asymptotic  approximation  to  I(k,8)  gives  good  accuracy,  and 
removes  the  need  for  numerical  tabulation. 

It  is  well-known,  and  was  shown  in  Appendix  B,  that  only  the  first 
Fresnel  zone  contributes  significantly  to  the  received  signal,  provided 
that  the  surface  contains  no  abrupt  nonuniformities  in  either  shape  or 
electrical  properties.  The  angular  width,  A$,  of  this  zone  is  of  order 

A*  * [A/Tr0l1/2  , (C-21) 

which  is  small  provided  that  the  transmission  pathlength  is  at  least  several 
wavelengths.  Thus,  subject  to  a posteriori  justification  we  use 

cos*=l-*2/2  , (C-22) 


in  Eq.  (C-18),  which  results  in  the  following  form  of  I(k,e): 


I = /27b 


eika/B(2(A-B)/B+*2]1/2 

t2(A-B)/B+*2]1/2 


(C-23) 


Further,  since  large  values  of  $ are  unimportant,  the  upper  limit  can  be 
changed  from  * to  •,  and  Eq.  (C-23)  can  be  approximated  by 


feu 


IsLin  /2/B 
b-M) 


fir- 


ika/¥[2(A-B)/B+*2]1/2 


? .■/, cosb$  , (C-24) 

I2(A-B)/B+<rr 


which  is  of  the  standard  form  (.Gradahteyn  and  Ryzhikt  1965) 


r ,ip[«V]1/2  nl 

J.  d*  “sb*  *'T 


hJ  (mVp2-b2  ) ■ (C-: 


where  the  positive  square  root  is  taken  and  Is  the  Hankel  function.  It 
follows  that 


I(k,0)  « ^ /2?B  hJ  (/2  ka  /a^b) 


(C-26) 


or,  after  using  Eqs.  (C-19)  and  (C-20), 


lm  * ¥ Vsine2sine.  M 2 ^ sin 


Vs 


1 a ‘T’  V r a2  , a hH  2 ka  sin 

2 TsinQ  sin0rt  0\ 


0 • 


(C-27) 


The  argument  of  the  Hankel  function  is  kr^  (see  Eq.  C-14) , evaluated  at 
$■0,  and  may  be  assumed  large  enough  to  use  the  asymptotic  expansion  of 
the  Hankel  function.  Thus,  to  order  1/ka  sin0,  substitution  of  Eq.  (C-27) 
into  Eq.  (C-17)  gives 
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W(s0) 


ira 


sin  Sg/2a 

sin 

i/a 

sin  s/2a  | 

^sin  Sg/a 

.JVl') 

V 2*  ) 

21ka(sln  s/2a-8in  Sg/2a  + |sin  (s-Sg)/2a 


D 


(C-28) 


where  a denotes  great-circle  distances  along  the  earth,  and  we  have  used 
8 ■ s/a,  etc.  For  s^Sq,  the  Integrand  oscillates  very  slowly,  as  can.  be 
seen  by  the  fact  that  the  exponent  vanishes  to  first  order  In  Sg/2a. 

For  s > 3q,  however,  the  Integrand  oscillates  very  rapidly,  as  can  be  seen 
from  the  fact  that  the  exponent  Is  2ik(s-Sg>  to  first  order  in  Sg/2a. 
Thus,  we  write  for  the  final  fora  of  the  Integral  equation  for  a smooth, 
round  earth 


W(s0)  * 


. m1/2  “iri/4  da_ 


W(s) 


sin  3g/2a 


sin  s/2a 


sin  s/a 


ksin  (Sg/a)  a ini 


2ika(sin  s/2a-sln  Sg/2a  + sln  (Sg-s)/2a) 


(C-29) 


which  is  Eq.  (11)  on  p.  11  of  the  main  text.  The  correction  terms  to  Eq. 
(C-19)  are  extremely  tedious  to  derive,  and  are  not  really  needed,  because 
the  excellent  agreement  between  numerical  solutions  of  Eq.  (C-29)  and 
available  results  from  the  residue  series  is  an  adequate  accuracy  check. 
Nonetheless,  we  point  out  that  by  following  steps  analogous  to  Eqs.  (B-13) 
through  (B-17),  we  find  after  considerable  algebra,  that  the  fractional 
error  of  Eq.  (C-29)  is  of  the  same  order  of  magnitude  as  the  error  term 
given  by  Eq.  (B-24) . 
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Appendix  D 

NUMURICAL  SOLUTION  OF  EQS.  (8)  AND  (II) 


HETHOD 

Both  Eq.  (8)  and  Eq.  (11)  (pp.  6 and  11)  are  Instances  of  the  Volterra 
integral  equation  of  the  second  kind: 

x 

W(x)  - g(x)  + J W(s)K(x,s)ds  . (D-l) 

0 

We  solve  this  aquation  iteratively  using  the  well-known  method  of  Picard. 
Starting  with  the  initial  guess  W^(x)  ■ g(x),  successive  approximations 
WL,  W2,  ...  are  found  such  that 


A 

Wj+i(x)  - g(x)  + y*  Wj(s)K(x,s)ds 


x fixed,  is  singular  at  the  endpoints  of  the  interval  (0,x) . Both  Eq. 
(8)  and  Eq.  (11)  can  be  rewritten  in  the  form* 


w(x)  ■ ,<*)  + f d. 


* 

This  procedure  is  possible  for  Eq.  (11)  because  sln(s/2a)  and 


sin  (~2^")  have  only  first-order  zeros  near  s«0  and  s*x. 


(D-2) 


Iteration  continues  until  two  successive  approximations  differ  by  less 
than  some  specified  tolerance.  Volterra  integral  equations  and  their 
solution  by  Picard's  iterative  method  are  discussed  by  Collatz  {I960). 

The  difficulty  with  the  solution  by  Picard's  method  lies  in  the 
repeated  evaluations  of  the  integral  in  Eq.  (D-l).  From  Eqs.  (8)  and 
(11),  it  is  seen  that  the  kernal,  K.(x,s),  viewed  as  a function  of  s with 


I 


(D-3) 


where  H(x,s)  has  derivatives  of  all  orders.  A Chebyshev-Gauss  quadrature 


is  used  to  evaluate  the  Integral  in  Eq.  (D-3) . For  each  iteration,  , 
ve  take 


f w.(s)H(x,s)  J1 


/s(x-s) 


(D-4) 


where 


fj(s)  - Wj  (s)H(x,  s) 


x x cos(2i-l)it 
12  2 2n 


(D-5) 


R f(2a>(5)  0 < 5 < x . (D-6) 

Q 2Zn(2n) ! 


Here,  Rq  is  the  remainder  term  corresponding  to  the  computational  error 
incurred  by  use  of  a finite  number,  n,  quadratures.  Note  that  this  remainder 
is  zero,  and  the  integration  exact,  when  f is  a polynomial  of  degree  less 
than  2n.  The  Chebsyhev-Gauss  quadrature  is  described  by  Krylov  (1962) . 

IMPLEMENTATION 

The  code  used  to  solve  Eq.  (8)  and  Eq.  (11)  is  written  in  FORTRAN. 

The  function  Wj , represented  as  a vector  array  of  chosen  length,  mathe- 
matically corresponds  to  representing  the  successive  approximations  to 
the  solution  as  a piece-wise  linear  function.  The  code  user  selects  the 
number  and  widths  of  the  steps  in  the  approximating  piece-wise  linear  func- 
tion. Iteration  according  to  Eq.  (D-2)  continues  until  a pair  and  W^+1 
are  found  that  satisfy  the  inequality 


( Max 

|E  [vv  - v(v] 


*r*m 
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where  N Is  the  number  of  values  in  the  vector  array  approximating  W, 

til 

x^  Is  the  value  of  x at  the  N step,  and  T is  a chosen  tolerance.  Thus, 

our  convergence  criterion  demands  that  the  rms  difference  between  two 

successive  iterations  be  less  than  the  chosen  tolerance. 

Most  of  the  results  given  in  the  main  text  are  obtained  with  a 

convergence  tolerance  of  10  , and  were  run  on  a CDC  7600  computer.  The 

running  time  required  to  calculate  W as  a function  of  distance  depends 

on  the  propagation  pathlength,  ground  conductivity,  wave  frequency, 

number  of  quadratures,  N , and  T.  The  running  time  tended  to  increase 

mflx 

as  the  frequency  increased  or  the  conductivity  decreased.  For  example, 
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for  T-10  , 11  sec  of  running  time  was  required  to  calculate  W out  to  a 
distance  of  2000  km  for  a conductivity  of  4 mhos/m  and  a frequency  of 
100  kHz.  Reduction  of  the  frequency  to  20  kHz  reduced  the  running  time 
to  8 sec;  but,  for  100  kHz,  decreasing  o to  2x10  ^ mhos/m  caused  18  sec 
of  running  time  to  be  required  to  calculate  W out  to  only  50  km.  Of  course, 
at  100  kHz,  about  as  much  attenuation  occurs  for  50  km  of  propagation  over 
ice  (2x10  ^ mhos/m)  as  for  more  than  1000  km  over  seawater  (4  mhos/m). 

Thus,  the  running  time  required  to  compute  a given  amount  of  attenuation 
does  not  vary  drastically. 


